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Abstract 

We consider the hydrogen quantum plasma in the Saha regime, where it almost reduces to a 
partially ionized atomic gas. We briefly review the construction of systematic expansions of thermo- 
dynamical functions beyond Saha theory, which describes an ideal mixture of ionized protons, ionized 
electrons and hydrogen atoms in their ground-state. Thanks to the existence of rigorous results, we 
first identify the simultaneous low-temperature and low-density limit in which Saha theory becomes 
asymptotically exact. Then, we argue that the screened cluster representation is well suited for cal- 
culating corrections, since that formalism accounts for all screening and recombination phenomena 
at work in a more tractable way than other many-body methods. We sketch the corresponding dia- 
grammatical analysis, which leads to an exact asymptotic expansion for the equation of state. That 
scaled low-temperature expansion improves the analytical knowledge of the phase diagram. It also 
provides reliable numerical values over a rather wide range of temperatures and densities, as confirmed 
by comparisons to quantum Monte Carlo data. 

1 Introduction 

Obtaining asymptotically exact formulae for the equation of state of quantum Coulomb matter 
is important, both at a theoretical level and for practical applications. They provide a better 
understanding of basic phenomena like molecular recombination and screening in the framework of 
statistical mechanics. Such formulae are free from any intermediate phenomenological modelization 
and uncontrolled approximation. They provide moreover reliable and accurate data in some range 
of thermodynamical parameters. Exact expansions are of particular interest for hydrogen described 
as a gas of quantum protons and electrons interacting via the Coulomb potential. Indeed, thanks 
to its relative simplicity, analytical calculations can be carried out further than for heavier species. 
In practice, the corresponding expansions are quite useful since hydrogen is the most abundant 
element in the universe. Astrophysicists need accurate equations of state over a wide range of 
temperatures and densities, including the so-called Saha regime where hydrogen reduces to a 
partially ionized atomic gas. 

In the 50 's pQ, a first kind of asymptotic expansion was derived for an electron gas at high 
densities, which behaves as a free Fermi gas in a first approximation. Corrections can be computed 
in a systematic way in the framework of standard many-body perturbation theory [2J, where the 
small-expansion parameter is the charge of the electrons. A second kind of asymptotic expansions, 
namely the familiar virial expansions, where constructed in the 60's [3], for two- or more component 
systems, including quantum hydrogen. In such expansions, temperature T is fixed at a non-zero 
value, and numerical densities p a for species a are driven to zero. For hydrogen, we set p — p p — p e . 
At lowest order, the system behaves as an ideal mixture of nuclei and electrons, as rigorously proved 
in Ref. [4] . Corrections are represented by series involving integer and half- integer powers of the 
Pq s, as well as logarithmic terms. Corresponding calculations have been first performed up to 
order p 2 a [5] , by using the effective-potential method introduced by Morita [B] . Further corrections 
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of order pd have been derived in the 90 's within another formalism based on the path integral 
representation [7], and retrieved later by Ebeling-Morita's method [8]. 
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Above expansions are suited for regimes where the systems are almost fully ionized. They 
cannot describe the Saha regime for hydrogen, where a finite fraction of electrons and protons 
recombine into atoms in their groundstate. Then, according to the familiar Saha theory [5], the 
system is expected to behave as an ideal mixture of ionized protons, ionized electrons and atoms 
H. The construction of suitable asymptotic expansions requires first to identify, if it exists, a 
regime of thermodynamical parameters where Saha theory becomes asymptotically exact. As 
described in Section 2, that regime is obtained by setting T — > 0, while ratio p/p* is kept fixed 
with temperature-dependent density p* given by expression 

Once the proper limit which defines the Saha regime has been identified, the construction of 
systematic expansions beyond Saha theory requires a formalism which accounts for all recombi- 
nation and screening phenomena at work. For that purpose, the screened cluster diagrammatical 
representation [10] is particularly adequate, as briefly described in Section 3. The corresponding 
analysis of all involved graphs, provides the so-called SLT expansion of pressure P, i.e. 



where Saha pressure Psaha is given by formula ([5]) in section 2. Functions bk(p/p*) only depend 
on ratio p/p*, while temperature dependent functions afc(/3) decay exponentially fast when T 
vanishes, ~ exp(— j36f.) except for possible multiplicative powers of (3. Expansion ([T]) is 

ordered with respect to the decaying rates < 8\ < 62 < — of the a/c(/3)'s functions. The first 
five corrections computed in Ref. [TT] are schematically presented in Section 4. They account for 
non-ideal phenomena such as plasma polarization, shift in the atomic energy levels, interactions 
between ionized charges and atoms, and also formation of molecules H2 or ions H~ and . 
Further corrections k > 6 decay exponentially faster than exp(/3i?^f), where Eh — — me 4 /(2h 2 ) is 
the atomic groundstate energy and m — m p m e / \m p + m e ) is the reduced mass for the two-body 
electron-proton problem. Along a given low-temperature isotherm, we also study the behaviour 
of the bk(p/ p*)'s. SLT expansion (fT]) then appears as a partial infinite resummation of ordinary 
virial expansions at low densities p <^ p*. Expansion (fT]) remains valid at intermediate (p ~ p*) 
and large (p > p*) densities, but breaks down at too large densities p ^> p* because molecular 
recombination becomes then prominent. 

As usual for asymptotic series, keeping only the first few terms of SLT expansion (TTJ) should lead 
to an accurate equation of state, provided that thermal energy fc#T is smaller than Rydberg energy 
\Eh\- In Section 5, we give a flavour of numerical calculations based on the truncation of (flj up to 
term k = 5 included [12]. Because of the relatively large temperature scale \Eu\/kB — 150000A", 
and of the occurrence of exponentially decaying factors, that truncated equation of state is reliable 
over a rather wide range of thermodynamic parameters, as confirmed by comparisons to quantum 
Monte Carlo simulations by Militzer and Ceperley [13] . This allows us to introduce a semi-empirical 
criterion for the convergence of SLT expansions, which provides the validity domain of the corrected 
EOS in the temperature-density plane. Our formulae should be particularly useful in physical 
situations where deviations to Saha theory play an important role even if they remain small. For 
instance, a very accurate EOS is needed for interpreting recent seismology measurements in the 
Sun [14]. Our corrected EOS should be useful since, according to usual models, the Sun adiabat 
lies in the previous validity domain (see Fig. [lj. However, notice that an accurate description 
of that adiabat requires to take into account heavier species like helium, carbon, nitrogen and 
oxygen. If exact calculations become much more complicated, a simple account of obvious ideal 
contributions, in particular for less abundant species, should be sufficient for significantly improving 
a pure- hydrogen EOS. 

We stress that our approach does not provide unambiguous definitions of neither free and 
bound charges, nor ionization rates. As commented in Section 2] this does not cause any trouble 
as far as thermodynamic quantities are concerned. Other quantities like conductivity or opacity 
cannot be obtained within the present formalism. Usually, such quantities are computed within the 
framework of the chemical picture, where atoms and molecules are introduced phenomenologically 
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as preformed entities. A first-principles description is of course possible in principle, but it becomes 
quite cumbersome since it would require a rigorous introduction of either real-time evolution or 
coupling to radiation. 



2 The hydrogen gas in the Saha regime 
2.1 The physical picture 

Within the physical picture, a hydrogen gas is viewed as a system of quantum point particles which 
are either protons or electrons, interacting via the instantaneous Coulomb potential v(r) = l/r. 
Protons and electrons have respective charges, masses, and spins, e p = e and e e = — e, m p and 
m e , cr p = a e = 1/2. In the present non-relativistic limit, the corresponding Hamiltonian for 
N = N p + N e particles reads 



N 



\ + -^2e ai e aj v(\xi-Xj\) (2) 
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where ai — p, e is the species of the ith particle and Aj is the Laplacian with respect to its position 
Xi. The system is enclosed in a box with volume A, in contact with a thermostat at temperature 
T and a reservoir of particles that fixes the chemical potentials equal to p p and p e for protons and 
electrons respectively. Because the infinite system maintains local neutrality p p = p e in any fluid 
phase, the bulk equilibrium quantities depend in fact solely on the mean 

p = {p p + Me)/2, (3) 

while the difference v = (p e — p p )/2 is not relevant as rigorously proved in ref. |15j . 

2.2 Identification of the scaled limit 

In the so-called Saha regime, a finite fraction of protons and electrons combine into hydrogen atoms, 
forming a partially ionized hydrogen gas. That regime is attained when several conditions are met. 
The temperature must be sufficiently low so that atoms can form, namely kT <C \Eh\- The density 
must be sufficiently low as well so that atoms maintain their individuality thanks to a 3> «b, 
where a is the mean inter-particle distance and cib is the Bohr radius. If the density becomes 
too low, atoms dissociate by entropy, while if it becomes too large, they recombine into molecules 
i?2- According to those simple considerations, both T and p must go to zero, in a related way. 
The precise form of that relation can be inferred from a rigorous analysis in the grand-canonical 
ensemble devised by Macris and Martin [TB] , who extended Fefferman's work on the atomic phase 
of the hydrogen plasma |17j . They introduce a scaling limit where the temperature T goes to 
zero, while the average chemical potential p of protons and electrons approaches the ground-state 
energy Eh with a definite slope [16] . Then, they proved that pressure P, within that scaling limit, 
tends to its Saha expression Psaha-, which describes an ideal mixture of hydrogen atoms, ionized 
protons and ionized electrons. 

In terms of temperature and density, the previous scaling limit can be rephrased as a low 
temperature expansion at fixed ratio p/p* [see eq. (1)], where p* is the temperature-dependent 
density 

= ™' h V = (W- (4) 

Notice that density vanishes exponentially fast when T is sent to zero. This ensures the proper 
energy-entropy balance which keeps a finite ionization rate that is entirely determined by the fixed 
ratio p/p*. Pressure P in units of p*ksT tends to Saha formula 

PPSaha/p* = P/P* + (1 + 2p/p*) 1/2 - 1 (5) 

apart from exponentially small terms when T — > 0. The respective behaviours of Saha pressure 
©, Psaha ~ 2/9 for p <^ p* , and Psaha ~ p for p 3> p* clearly illustrate that p* is a cross-over 
density between full ionization and full recombination. 
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3 Construction of SLT expansion 



3.1 Introduction of a suitable formalism 

Corrections to Saha pressure ([5]) involve interactions between ionized charges and atoms, as well 
as formation of ions and molecules. Standard many-body theory is not well-suited for taking into 
account recombination, since that mechanism is not pcrturbative with respect to the charge. For 
instance, an infinite number of Feynman ladder graphs must be resummed for describing a single 
atom H . 

The effective-potential method, which amounts to introduce a classical equivalent system of 
point particles with many-body effective interactions, is a priori more efficient for dealing with 
recombination. Indeed, n-body effective interactions are inferred from n-body quantum Gibbs 
factors which do account for recombination of n particles at short distances. The contributions 
of two-body effective interactions can be analyzed within standard methods of classical statis- 
tical mechanics. In particular, that feature has been exploited for computing virial expansions 
up to order p 2 [3], where atomic contributions appear. Ionic and molecular recombination are 
embedded in three- and four-body effective interactions. Unfortunately, the analysis of the cor- 
responding contributions becomes rather cumbersome, in particular because no standard classical 
tool is available. 

Above drawbacks of both standard many-body theory and effective-potential method, clearly 
emphasize the need for a formalism more appropriate to deal with Saha regime. The ACTEX 
method, introduced by Rogers [TH] in the 70's, is intended to account for the formation of chemical 
species in the framework of the physical picture. That approach starts from the usual activity 
expansion of thermodynamical quantities in the grand-canonical ensemble. Despite its rather 
successful predictions at moderate densities and temperatures, that approach cannot be applied 
here as it stands, because ACTEX series are not exactly reorganized via a systematic treatment of 
both recombination and screening. In fact, quantum Gibbs factors involved in activity series do not 
factorize as products of two-body counterparts, like in the case of classical charges. Resummations, 
which are crucial for taking into account screening effects, are very hard to handle in the usual 
quantum activity series. 

The difficult task of controlling simultaneously recombination and screening effects in quan- 
tum activity series can be accomplished thanks to the Feynman-Kac path integral representation. 
Within that formalism, the genuine quantum system of point protons and electrons is shown to 
be equivalent to a classical gas of extended loops [19]. Thermodynamic quantities of hydrogen 
are then represented by activity series in the world of loops, which can be suitably rearranged as 
described below. 

3.2 The screened cluster representation 

Since loops are classical objects with two-body interactions, standard Mayer diagrammatical meth- 
ods can be applied. In particular, the one-body loop density is represented by a series of Mayer-like 
graphs in the grand-canonical ensemble. Usual points are replaced by loops, loop fugacities are 
simply related to particle fugacity z = exp(/3/z), while Mayer bonds are built with the loop- loop 
interaction potential. Since that potential behaves as the Coulomb interaction at large distances, 
Mayer graphs are plagued with long-range divergences. Such divergences are systematically re- 
moved via chain resummations, which amount to introduce a screened potential [20] • Contrary 
to the familiar classical Debye potential which decays exponentially fast, that quantum potential 
decays only as 1/r 3 at large distances r. However, at low densities, it reduces to its classical Debye 
counterpart plus small corrections. At the same time, the whole series is reorganized in terms 
of particle clusters. Eventually, particle density p, obtained by integrating loop density over all 
possible shapes, is exactly rewritten as the following diagrammatical series [10] 
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The corresponding graphs are constructed with topological rules close to that of ordinary Mayer 
graphs, except for some exclusion constraints avoiding double counting. Usual points are now 
replaced by particle clusters. The statistical weight of a given cluster involves particle fugacities, 
as well as screened interactions. Two clusters can be connected by a single bond, which is either 
— /3<5, /3 2 $ 2 /2, or — /3 3 $ 3 /6, where 3> is the screened interaction between those clusters. 

Expansion ([6]) accounts in a fully consistent way for all effects of interactions in the system 
at finite density and finite temperature. The various phenomena at work are embedded in well- 
defined graphs. For instance, the first graphs shown in © account respectively for a single ionized 
proton, formation of an atom H and of a molecule H2, and interactions between two atoms. As 
the scaling limit defines quite diluted conditions, only a few simple graphs in the screened cluster 
expansion are expected to contribute to the first corrections to Saha theory. 

3.3 Behaviour of graphs in the scaling limit 

The scaling limit defined in Section 2.2 within the grand-canonical ensemble, can be rephrased as 
follows. Starting from particle fugacity z and temperature T, we introduce a new couple (7, T) 
of independent thermodynamic parameters defined through the relation z equal a constant times 
7 exp(/3.Eff). Then, we set T — > at fixed 7. The behaviours of graphs in representation ^ result 
from the competition between three mechanisms, which can be roughly described as follows. 

• Screening : Contributions of bonds — /3<£>, /3 2 $ 2 /2 and — /3 3 $ 3 /6 are controlled by the inverse 
Debye screening length k = (87r/3e 2 p*7) 1 / 2 for ionized protons and ionized electrons with density 
p*7. They behave as positive or negative powers of n, which itself decays exponentially fast as 



• Recombination : For a particle-cluster made with N p protons and N e electrons, the behaviour 
of its statistical weight gives raise to a cluster partition function Z(N p ,N e ) in the vacuum. Each 
Z(N p ,N e ) is a truncated trace of Gibbs operator exp(— /3_ffjv p .ArJ involving only bare Coulomb 
Hamiltonians, which converges thanks to a truncation inherited from screening by ionized charges |10) 
In Z(N p , N e ), contributions of all possible recombined entities made with M p < N p protons 
and M e < N e electrons, are mixed together. Remarkably, the contribution of a given chemi- 
cal species made with N p protons and N e electrons, naturally emerges through Boltzmann factor 
exp(-/3£W ), where £$> N is the groundstate energy of Hamiltonian Hn p ,n c - That factor 
increases exponentially fast since E$ N < 0. 

• Entropy : In a given cluster, the presence of N — N p + N e particles generates activity powers 
z , which decay exponentially fast as exp(Nf3En). 

The behaviour of a graph in the scaling limit is obtained, roughly speaking, by taking the 
product of the exponential factors generated by each of the above mechanisms. Then, every graph 
is found to decay exponentially fast. The leading contributions arise from the first two graphs 
in representation ([6]). They do reduce to the ideal terms p l p and p 1 ^ predicted by Saha theory. 
Further corrections decay exponentially faster than p* in agreement with rigorous bounds [15] . 
Dividing all terms by p*, the SLT expansion of p/p* reads [TT] 



where the first two terms account for contributions from ionized particles and hydrogen atoms in 
their ground state. The remaining terms involve functions /ifc(/3) that decay exponentially fast 
when T — y 0, while r f lk is an integer or half-integer power of 7 (which may be multiplied by integer 
powers of In 7 when k > 6). Expansion ([7]) is ordered with respect to increasing decay rates of 
the hkS. The corresponding hierarchy follows from subtle inequalities between the groundstate 
energies of all Coulomb Hamiltonians Hn v ,n c - For instance, the molecular contribution, which 
determines the leading low-temperature behaviour of /12, indeed decays exponentially fast thanks 
to Eh 2 > 3Eh '■ this ensures that molecules H2 are very scarce in the Saha regime compared to 
atoms H, despite they are more stable energetically, i.e. Eh 2 < 2Eh. 



exp((3E H /2). 




(7) 
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4 Equation of state beyond Saha theory 



4.1 Scaled low-temperature expansion of pressure 

Representation ([7J expresses the density in terms of variables T and 7, or equivalently T and 
p since there is a one-to-one correspondence between those sets of variables. As the natural 
thermodynamical parameters are the temperature and the density, it is quite useful to invert the 
SLT expansion ([Jj to determine j(p,T), in order to compute all thermodynamical quantities as 
functions of T and p. In the present scaling limit, this can be done in a perturbative way around 
Saha expression 

lSaha(p, T) = VT+W - 1, (8) 

easily obtained by keeping only the first two terms in ([7]). Each correction to that form reduces to 
a product of an algebraic function of pj p* times a temperature-dependent function which decays 
exponentially fast. 

The standard thermodynamical relation, which expresses density in terms of the partial deriva- 
tive of pressure with respect to z at fixed T, can be rewritten here as 

0.) 

After inserting SLT expansion ([7]) of p into ©, a straightforward integration with respect to 7 
provides the SLT expansion of P in terms of 7 and T. The corresponding expansion (fTJ) of f3P/ p* 
in terms of pjp* and T, then follows by using the inversion relation y(p,T) determined above. 
The physical content of first five corrections in (JTJ), as well as the expressions and values of the 
corresponding decay rates are summarized in the following table 



Correction (k) 


Physical content 


S k (in eV) 


1 

2 
3 
4 
5 


plasma polarization around ionized charges 
formation of molecules, atom-atom interactions 
atomic excitations, charge-charge interactions 
formation of ions, atom-charge interactions 
fluctuations of plasma polarization 


|£ff|/2~6.8 
\3E H -E H2 \~9.1 
3|£ H |/4~ 10.2 
\2E H -E H +\ ~ 11.0 
\Eh\ ~ 13.6 



First correction k = 1 is equivalent to a modification of Saha ionization equilibrium |21) derived 
within Green functions techniques (see also Ref. [3]), where rate 81 = \Ejj\/2 arises from the 
behaviour k ~ exp(— /3\Eh\/2) in the scaling limit. All further corrections are entirely new, as 
well as the structure of SLT expansion (p}. Beyond their leading behaviours «fe(/3) ~ exp(— /38k), 
functions aik{/3) include further corrections which decay exponentially faster. For instance, if the 
leading behaviour of 013 (/3) is controlled by the first atomic excited state, all the other contributions 
of excited states are incorporated into a^{j3). Similarly, «2(/3) includes not only the contribution 
of the molecular groundstate, but also all contributions of molecular excited states. As mentioned 
above about recombination, the sums of all those contributions are indeed finite. Also, we stress 
that such contributions of recombined entities, like atoms H in k — 3 or molecules H2 in k = 2, are 
entangled with that of their dissociation products. Thus, purely atomic or molecular contributions 
cannot be unambiguously defined. This does not any trouble here, since only full contributions 
embedded in k = 3 and k = 2 are relevant for thermodynamics. In other approaches based on the 
chemical picture, that ambiguity has been the source of many controversies since the introduction 
of Planck-Larkin formula (see e.g. Ref. [22] . Ref. [23] and Ref. [24]). Eventually, notice that the 
contributions to expansion ([1]) of more complex entitites, like , or H3, decay exponentially 
faster than exp(— (3\Eh\) as detailed in Ref. [TTj . 

4.2 Low-temperature isotherms 

Let us consider now a small fixed temperature T, and study the behaviour of various corrections to 
Saha pressure in (fT]) when density p is varied. The corresponding low (p <§; p*) and large (p 3> p* 
) density behaviours are summarized below : 
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k 


p<- p* 


P> p* 


1 


p 3/2 


p 3/4 


2 


P 4 


P 2 


3 


P 2 


p l/2 


4 


P 3 


p 3/2 


5 


P 2 


p l/2 



At low densities p -C p*, the leading correction of order p 3 / 2 is given by plasma polarization 
(term k = 1), while at large densities p ^S> p* the leading correction of order p 2 arises from both 
molecules i?2 and atom-atom interactions (term k = 2). We have checked that the familiar virial 
expansion at low densities is indeed recovered up to order p 2 included, since all terms k > 6 provide 
powers higher than p 2 . We stress that at too large densities, expansion (TTJ breaks down because 
various corrections, in particular those due to molecular recombination, prevail over Saha pressure 
which grows only as p. In fact, a semi-empirical criterion based on that observation allows us to 
infer a validity domain for the SLT expansion, as described below. 

5 Numerical applications and comparisons to Monte Carlo 
data 

The truncated EOS obtained by keeping the first five terms in SLT expansion (|TJ) can be computed 
numerically. The b^s are easily computed since they reduce to simple algebraic functions of 
ratio p/p*. Numerical values for temperature-dependent functions a±, a^, and a§ are also readily 
inferred from explicit analytic expressions. No similar expressions for a-i and are available, 
since analytical results on the three- and four-body quantum problem are very scarce. Then, we 
use simple modelizations of those functions which account for their exact low-temperature forms 
on one hand, and incorporate familiar phenomenological descriptions of ions H~ and H% and of 
molecule H2 on the other hand [12] . 

Various isotherms corresponding to increasing temperatures have been considered. For tem- 
peratures below 2000A, the Saha regime defines extremely diluted conditions which do not make 
physical sense : this explains why Earth or Brown Dwarfs atmospheres only involves molecules Hi- 
Temperature can be increased up to T = 30000A, which is still small compared to the characteristic 
temperature scale \Eff\/kg ~ 150000A. At a given temperature, calculations within the truncated 
EOS are not reliable above some density p c , for which corrections to Saha pressure become too 
large, in agreement with previous estimations for p 3> p* ■ That breakdown is due to molecular 
recombination below T ~ 16000 K, and to atom-atom interactions for higher temperatures. 

At relatively low temperatures, i.e. below T = 10000A, PIMC calculations [13] have been per- 
formed at rather high densities p » p c for which atoms H are mainly recombined into molecules 
H2- In the Saha regime, statistics in PIMC results are poor because the corresponding densities 
are too diluted. In fact, under such conditions, our analytical results might serve as testbenchs 
for simulation methods. For higher temperatures, i.e. above T — 10000A", there exists a density 
range (see PIMC crosses displayed in Fig. validitydomain) , where comparisons between our re- 
sults and PIMC data [13] are instructive. A good agreement is observed in some density range, 
which can be inferred from a semi-empirical criterion : corrections cannot exceed a few per cent 
of Saha pressure. This defines, at a heuristic level, the validity domain of SLT expansion shown 
in Fig. validitydomain. In that whole domain, weak-coupling and weak-degeneracy conditions 
are fullfilled. The tongue structure of the domain between T ~ 10000 K and T ~ 25000 K results 
from the increase of the strength of interactions between ionized charges and atoms. Notice that 
the whole domain is restricted to rather low densities in general, so high-density phenomena, like 
the celebrated plasma phase transition, remain beyond the scope of our approach. 

Eventually, we emphasize that our numerical calculations will be detailed in a forthcoming 
paper [12j . where simple representations of functions ak(/3) (k = 1, 5) will be given, while the 
corresponding functions bk(p/p*) can already be found in Ref. [11]. Our results will be compared 
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to PIMC data, and also to phenomenological calculations. Applications to the Sun adiabat should 
be considered later. 
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